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IMPROVED EVENT LOCATION UNCERTAINTY ESTIMATES 

Istvan Bondar 1 , Keith McLaughlin 1 , and Hans Israelsson 1 
Science Applications International Corporation 
Sponsored by Air Force Research Laboratory 
Contract No. FA8718-05-C-0018 


ABSTRACT 


While many recent studies aimed to reduce location bias by introducing improved travel-time corrections, less effort 
was devoted to the complete estimation of location uncertainty, despite the fact that formal error ellipses are often 
overly optimistic. Since most location algorithms assume that the observations are independent, correlated 
systematic errors that are due to similar ray paths inevitably result in underestimated location uncertainties. 
Furthermore, the tails of real seismic data distributions are heavier than Gaussian. The main objectives of this 
project are to develop, test, and validate methodologies to estimate location uncertainties in the presence of 
correlated, systematic, and non-Gaussian errors. Particular attention is paid to robust and transportable models for a 
travel-time covariance matrix. 

To address correlated errors, we estimate the spatial correlation structure in arrival-time data using variogram 
models. We developed a methodology based on copula theory to derive robust, data-driven variogram models. For 
validation purposes, we use GTO-2 event clusters. These include the Nevada, Lop Nor, Semipalatinsk, and Novaya 
Zemlya test sites, as well as the Azgir Peaceful Nuclear Explosions and the Lubin, Poland, mine-related events. 
Using ground-truth (GT) clusters allows us to calculate “ground truth” residuals with respect to the GT locations for 
a specific velocity model. We show the improvements in variogram estimates when using global 3D models instead 
of the iasp91 model. The proper choice of the underlying velocity model is especially important for regional phases. 

To address issues embodied by real data distributions, we performed fully controlled experiments using known high 
signal-to-noise ratio (SNR) waveforms, scaled down to several magnitude levels and embedded in clean noise to 
derive models of measurement errors. We model these measurement errors as a series of generalized extreme value 
distributions whose parameters (location, scale, and shape) depend on the measured SNR. This allows us to model 
the increasing picking error bias and the increasingly heavier tails of the residual distributions with decreasing SNR. 

We have incorporated the full covariance matrix estimate in a linearized location algorithm. We show that by taking 
into account the correlated error structure with a robust transportable station-station correlation model, w r e achieve 
90% coverage (i.e., the 90% error ellipse covers the true location 90% of the time), for sparse, unbalanced networks. 
Furthermore, the coverage statistics do not deteriorate with an increasing number of stations. Future work will 
incorporate SNR-dependent, non-Gaussian arrival-time measurement errors. 
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OBJECTIVES 


The objectives of this project are to develop methodologies to estimate location uncertainties in the presence of 
correlated, systematic model errors and to characterize measurement errors as a function of signal parameters such 
as phase and SNR. The improved understanding of the complete error budget described by the full covariance 
matrix is incorporated into a linearized location algorithm, leading to more robust estimates of location uncertainty. 
The ultimate goal of this project is to develop transportable error models that will provide reliable location 
uncertainty estimates for small events recorded only by a few stations. 

RESEARCH ACCOMPLISHED 

The motivation for this project is the observation that location uncertainty estimates are often underestimated, that 
is, the error ellipses scaled to the 90% confidence level do not contain 90% of the true locations. While most recent 
location calibration studies focused on producing improved travel-time predictions to reduce location bias, less 
effort was devoted to obtaining reliable formal error ellipses. Linearized location algorithms typically rely on the 
assumption that the observations are independent and the residuals are Gaussian distributed. Unfortunately, these 
two assumptions are almost always violated. Many researchers have pointed out (e.g., Buland, 1986; Anderson, 
1982) that the distribution of measurement errors is non-Gaussian and better described by non-zero mean, skewed, 
heavy-tailed distributions. Furthermore, station distributions are far from uniform, and observations are not 
independent. Closely spaced stations sample similar ray paths and thus introduce correlated systematic errors. This 
redundancy in the observ ations reduces the effective number of degrees of freedom, thus ignoring the correlated 
structure inevitably results in underestimated location uncertainty estimates. 

In this project we focus on the treatment of correlated errors combined with non-Gaussian, non-zero-mean, heavy¬ 
tailed, skewed distributions of reading errors. Figure 1 illustrates our research strategy. Both the measurement and 
model errors are described by their covariance matrices, where the full data covariance matrix ( C D ) is represented by 
the sum of the reading error ( C R ) and network (C/v) covariance matrices. The network covariance matrix describes 
the correlated error structure and is estimated by variogram analyses. The reading errors are modeled by a general 
extreme value distribution, with parameters slowly increasing with decreasing SNR. The full data covariance matrix 
is incorporated in a linearized inversion scheme. Because of the non-Gaussian error distributions, the linearized 
inversion provides only an approximation of the solution; we will develop a hypothesis test to justify the validity of 
the “linearized” estimates. 
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Figure 1. Progress and strategy of research. Gray indicates the current state of the art linearized location 
algorithm. Validation test results are presented in this paper for the green stage (recent progress); 
yellow and orange represent planned future development. 

All the methodologies developed in the course of this project are tested and validated using GT event clusters, 
including GT0-2 underground nuclear explosions from Yucca Flat and Pahute Mesa, Balapan and Degelcn 
mountains, Novaya Zemlya, Azgir, and Lop Nor, as well as mining explosions from the Lubin mines. 
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Network Covariance Matrix 

The network covariance matrix accounts for the spatially correlated travel-time structure that is due to similar ray 
paths. We use variogram analysis to estimate the correlation structure in the data. In order to obtain robust 
variogram estimates, we use an entire GT event cluster. The variogram in geostatistical analysis is defined by the 
equation 

Y (h) = ^(&(A) - <5t( A + h)) : ^ = ct* u - Corr(h)CT^ = ct^ - Cov(h), (1) 

where o] m denotes the background variance, <5/(A) and dt(A+h) are the GT residuals (residuals with respect to 
the GT location, using a specific Earth model) at station pairs of common events, and h is the station separation. We 
employ copula formalism (Nelsen, 1999; Frees and Valdez, 1998; Genest and Rivest, 1993) to determine robust 
variogram models. A copula couples marginal distributions to form the joint probability distribution: H(x x v) = 
C(F(jcT G(v)) In other words, the copula describes the joint distribution of the order statistics u = F(h) and 
v = G M <*X(A) - At(A + h)J LWe have evaluated a number of copulas and found that the Clayton copula, defined 
as C(u}v) = (u" a + v' a - 1 j , a E(0,<»), provides the best fit to the data (Bondar et al, 2005). We define the robust 
variogram model as the median regression curv e of the residual difference squares for station pairs of common 
events with respect to station separation. Using the copula formalism, the median regression curve simply becomes 
the solution of the equation 

C'( v 111) = dC( d "’ V> - 0 5 , which yields y (h) -G' 1 ([l + iT°(o.5* a/(at,) - l)j’‘ “\,u £[0,1]. (2) 

Thus, the copula formalism offers a fully data-driven approach to derive robust, monotonically increasing variogram 
models that are free from Gaussian assumptions and provide analytic expressions parameterized by the order 
statistics. Having derived a variogram model, the estimates for the elements of the network covariance matrix are 

C N ( i ’j) = <^ 1| -y(A(sta j ,sta.)) . (5) 

Note that the isotropic variogram model defined by Eq. (1) does not account for azimuthal variations in the 
correlation structure. Hence, the proper choice of the underlying velocity model to calculate residuals becomes 
important, especially for regional phases. Figure 2 shows the Pn variograms obtained for Yucca Flat, Azgir, and Lop 
Nor GT event clusters when using iasp91 (Kennett and Engdahl, 1991) predictions to calculate GT residuals. The 
variance of the residuals, indicated by the binned median (blue line), increases significantly at far regional distances, 
where stations from different tectonic provinces are mixed together. This is poorly modeled by the isotropic 
variogram model, shown by the red line. 



Figure 2. Pn variograms for three GT event clusters using iasp91 predictions. The blue lines represent the 
median of every 5 percentiles; the red lines indicate the copula variogram models. Because of the 
unmodeled velocity structure by the iasp91 model, the isotropic variogram model provides a poor 
approximation of the correlation structure. 


To remedy these shortcomings, we recalculate the GT residuals using the global upper-mantle CUB2 model 
(Shapiro and Ritzwoller, 2004). The CUB2 model was extensively validated in Eurasia (Ritzwollcr et al., 2003; 
Yang et al., 2004), but here we show an example for Pn observations from Yucca Flat events. Figure 3 shows the Pn 
travel-time correction surface relative to iasp91 predictions and the comparison of the iasp91 and CUB2 GT 
residuals as a function of epicentral distance. The CUB2 residuals are closer to zero mean than the iasp91 ones, 
indicating that the CUB2 significantly reduces the path effects that were left unmodeled by iasp91. 
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Figure 3. (a) CUB2 Pn travel-time correction surface to iasp91 centered on Yucca Flat (star). The 289 stations 
with Pn readings are shown as triangles, (b) GT residuals using iasp91 (blue) and CUB2 (red) as a 
function of epicentral distance. 


Figure 4 shows the Pn variograms for the same event clusters as above, but now using the CUB2 predictions. The 
isotropic variogram models not only fit the observations better but also reflect a 50%-60% variance reduction 
typically achieved by the CUB2 model over iasp91 



Figure 4. Pn variograms for three GT event clusters using CUB2. The blue lines represent the median of 
every 5 percentiles, the red lines indicate the copula variogram models. Using a calibrated Earth 
model, the isotropic variogram models provide acceptable descriptions of the correlation structures. 


Figure 5 shows the comparison of the Pn variogram models derived for the various GT event clusters when using 
iasp91 or CUB2 travel-time predictions. The variance reduction provided by the calibrated CUB2 predictions is 
reflected in the reduced and remarkably consistent background variance (sill). 




Figure 5. Pn variogram models for GT event clusters using (a) iasp91 and (b) CUB2 predictions. The CUB2 
model produces much more consistent variograms, with reduced background variance (sill). 
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Modified Location Algorithm 

Standard linearized location algorithms assume independent Gaussian errors and solve the inversion problem by an 
iterative, weighted least-squares algorithm by minimizing the expression (</ -GrnJC R (c/ -G///), which is 
equivalent to solving the equation WGm = Wd, where G is the ( NxAf) design matrix containing the travel time 
derivatives for an event-station path, m is the ( Mxl ) model adjustment vector [a t , Ay, A v, AzJ , d is the ( Nxl) 

vector of time residuals, and W = C R 2 is the diagonal ( NxN) weight matrix. The G w m = d w equation is often 
solved by singular value decomposition, which yields G # = V W A~^U^. and thus the model adjustment of 
m est m G#d w . At each iteration the model vector is adjusted such that m j ^ 1 = w 7 + m^,. Once a convergent solution 
is obtained, the location uncertainty is defined by the a posteriori model covariance matrix, 

C M =G~ l C R G~ lT = V W A , which is ty pically scaled to the 90% confidence level. This linearized location 
algorithm constitutes our baseline, against which we measure improvements in location uncertainty estimates. 

In the presence of correlated systematic errors, the data covariance matrix is no longer diagonal. Our modified 
linearized location algorithm seeks a transformed set of equations WGm = Wd , in which the data covariance matrix 
is diagonal. To ensure this, we solve the inversion problem in the eigen coordinate system in which the transformed 
observations are independent. The singular value decomposition of the full data covariance matrix is written 

as C D = U D A D Vl , where A D is the diagonal matrix of eigenvalues and the columns of U D contain the eigenvectors 
of C D . We keep the first p largest eigenvalues from the cumulative eigenvalue spectrum such that 95% of the total 
variance is explained: ^A y /^A, £ 0.95 . We then define p as the effective number of degrees of freedom of the 
data, with an N-p dimension null space. Note that the 95% total variance level is somewhat arbitrary but a 
conservative and workable choice. Let C D = BB T , with B = U p A ll p 2 , then the projection matrix W = B - i =a -; /2 g; 
orthogonalizes the data set and projects redundant observations into the null space. After applying the projections 
G w * A~p 2 U T p G and d w = A~ p l2 U T p d , the formalism remains the same as in the baseline algorithm, but now d» 
represents linear combinations of the observ ed residuals, the “eigen residuals.” 



Figure 6. Upper panel: Network covariance matrix for a Lubin rockburst and its corresponding eigenvalue 
spectrum (red line). The green line shows the cumulative sum of eigenvalues as a percentage of the 
total sum. Only 8 eigenvalues are needed to explain 95% of the total variance. Bottom panel: Full 
data covariance matrix and its eigenvalue spectrum. Random picking errors blur the correlation 
structure and reduce the data redundancy. 
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The process is illustrated for a Lubin, Poland, event that occurred on May 26, 1995. Pn was reported at 96 stations, 
representing a dense but rather unbalanced network, as shown in the Figure 6 inset. The upper panel of Figure 6 
shows the network covariance matrix and its corresponding cumulative eigenvalue spectrum. The network 
covariance matrix has been arranged by its nearest-neighbor ordering of stations, and it exhibits a quasi-block- 
diagonal structure. Because many observ ations are strongly correlated, the first 8 largest eigenvalues explain 95% of 
the total covariance. In other words, 92% of the data carry redundant information. When the measurement (picking) 
error covariance matrix (assuming uniform, 0.5 s reading errors) is added to the network covariance matrix (lower 
panel), the reading errors weaken the correlation strength but leave the correlation pattern unchanged. Because of 
this “blurred" correlation structure, more eigenvalues are needed to explain 95% of the cumulative variance. 

In general, when the picks suffer from large measurement errors, the 
correlation structure is blurred by the random noise of picking errors. 

However, when the onsets are picked very accurately, the correlation 
structure becomes very important for obtaining reliable location 
uncertainty estimates. Figure 7 shows location results with the 
assumption of independent observations compared with those 
accounting for correlated errors. When the correlation structure is 
accounted for, not only does the error ellipse cover the GT location 
but the mislocation is also reduced from 15.6 km to 5.5 km. 


Validation Tests 

To validate the improved location uncertainty estimates that are due 
to the full data covariance matrix, we carried out two experiments 
using a set of events from each GT event cluster. In each experiment 
we located the events by both the baseline (independent error 
assumption) and the modified (correlated errors) location algorithms. 

We assume uniform, 0.5 s standard deviation measurement errors. For errors (red). 

Pn we use the CUB2 travel-time model and construct the network covariance matrices from the variogram models 
obtained with CUB2 predictions. For teleseismic P we use iasp91 travel-time predictions. 

The first experiment is designed to measure the robustness of the location uncertainty estimates against increasing 
levels of data redundancy. We locate events with an increasing number of stations of optimal azimuthal coverage 
and measure the mislocation, the area of the 90% error ellipse, and the GT coverage (whether the error ellipse 
contains the GT location). Figure 8 summarizes these results. When the correlated errors are accounted for, the 
improved relative weighting scheme reduces location errors. Since the inversion is performed in the eigensystem, 
redundant data are projected out and can no longer conspire to increase mislocation. The remaining mislocation 
represents the bias introduced by imperfect travel-time predictions. Furthermore, as the information content carried 
by the network is exhausted, the size of the 90% error ellipse stabilizes, while maintaining GT coverage. Thus, 
incorporating the correlated data structure in the location algorithm provides robust location uncertainty' estimates 
unaffected by an increasing degree of redundant information. 



uncertainties of a Lubin GT1 
event (star) using independence 
assumption (blue) and 
accounting for correlated 





Figure 8. (a) Mislocation, (b) error ellipse size, and (c) GT coverage with an increasing number of stations for 
an optimal network. Baseline location algorithm (independent error assumption) results are shown in 
blue; modified location algorithm (correlated errors) results are shown in red. 
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The second experiment is designed to measure the trustworthiness of the location uncertainty estimates for 
suboptimal networks. We locate events with randomly selected networks of 5, 8, 10, 20, and 50 stations (1000 
realizations each) and measure the distribution of the coverage parameter (normalized distance inside vs outside the 
90% ellipse). Note that under Gaussian assumptions, the coverage parameter follows a x 2 distribution with 2 degrees 
of freedom. The actual GT coverage is the percentile value, where the coverage parameter equals to one. For honest 
90% confidence error ellipses, this should obviously occur at the 90th percentile. Figure 9 shows the results for 
sparse and dense networks. While the assumption of independent observations produces abysmal coverage, taking 
into account the correlated error structure maintains 90% coverage for sparse, unbalanced networks. For dense 
networks, the coverage falls somewhat below 90%, and the deviation from the theoretical x distribution may 
indicate residual unmodeled nonlinear dependence structures or non-Gaussian effects. 



distribution of coverage parameters is show n in blue w hen events are located assuming independent 
errors and in red when using the full data covariance matrix. The theoretical distribution is drawn in 
black. When correlated errors are accounted for (red), the 90% confidence ellipse covers the true 
location for 90% of the realizations. 

Preliminary Model of Measurement Errors 

Measurement errors are typically modeled as Gaussian, zero-mean processes. However, picking errors of onset 
times suffer from heavy tails (Buland, 1986) and are skewed, as weak emergent arrivals are typically picked too late 
(Anderson, 1982). Measurement errors also suffer from systematic errors, as onset times along the same ray paths 
are systematically picked late with decreasing event size, or more precisely, with decreasing SNR (e g., Douglas et 
al., 1997, 2005). Douglas et al. (2005) point out that automatic detections are more likely affected by the systematic 
errors than are manual picks made by experienced analysts. Systematic reading errors introduce location bias, 
especially for small events recorded by sparse, unbalanced regional networks. 

To develop models of measurement error bias and variance, we follow the methodology developed by Kohl et al 
(2004, 2005). The methodology scales known, high-SNR signals (explosions and/or earthquakes) down to various 
magnitude levels and embeds them in clean background noise, thus enabling large-scale controlled experiments. We 
generated some 22,000 realizations of earthquake waveforms from the Lop Nor region scaled to the range of mb 
2.5-6.0, as well as some 9,000 realizations of Lop Nor explosion waveforms (mb 2.5-5.7), each embedded in 
different realizations of clean noise. Since the onset time and the SNR of the embedded signals are exactly known, 
running a signal detector allows us to derive improved measurement error models conditional on SNR. Figure 10 
shows the delays of the automatic picks w ith respect to the true onset times for first arriving phases with SNRs 
larger than 3.5. It is apparent that the onsets are picked increasingly late, with decreasing SNR. Furthermore, the bias 
and the scatter around the bias are much smaller for the more-impulsive explosion signals, and picking errors on the 
earthquake signals exhibit heavier tails. 
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a) 


i, 
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SNR 



Figure 10. Delay of automatic first-arriving P picks relative to true onset times for (a) earthquake and (b) 

explosion waveforms. Picking errors of earthquake seismograms exhibit larger bias and scatter and 



suffer from heavier tails than those from explosions. 



To model this increasing bias, scatter and 
skewness with decreasing SNR, we fit a 
series of general extreme value (GEV) 
distributions to the residual distribution 
as a function of SNR. Figure 11 
illustrates the preliminary, SNR- 
dependent measurement error models for 
the earthquake and explosion picks. The 
location, scale, and shape parameters of 
the GEV distributions vary slowly with 
SNR, allowing for continuous change in 
reading errors. The increasing reading 
error bias with decreasing SNR (dashed 
lines in Figure 11) introduces systematic, 
correlated picking errors for a fixed 
event. However, once the bias is 
corrected for, the reading errors become 
nearly independent; thus, the reading 
error covariance matrix may still be approximated by a diagonal matrix. 


Figure 11. Measurement error models for (a) earthquake and (b) 
explosion signals for first-arriving P presented as a series 
of GEV distributions, w ith parameters slow ly changing 
with SNR. 


CONCLUSIONS AND RECOMMENDATIONS 


We have developed methodologies to obtain robust estimates of network covariance matrices through copula-based 
variogram models. We have shown that the CUB2 global upper-mantle model provides significant variance 
reduction over the iasp91 model for regional Pn and thus produces reliable isotropic variogram models. We have 
developed a modified linearized location algorithm to incorporate the full, non-diagonal data covariance matrix in 
the inversion problem. Validation tests demonstrated that when the correlated data structure is taken into account, 
we 


• obtain robust location uncertainty estimates, unaffected by the amount of redundant information, 

• reduce mislocation that is due to conspiring errors, and 

• achieve 90% coverage 90% of the time for both sparse and dense unbalanced networks. 

While the variogram models for the various GT event clusters produce slightly different covariance matrices, the 
overall correlation structures are similar as the correlation varies the most at relatively small station separations (up 
to a few hundred kilometers). Preliminary tests indicate that it is possible to develop generic transportable variogram 
models for both Pn and P that will err on the conservative side and may work well anywhere on the globe. 
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We have developed preliminary measurement error models for first-arriving phases. We model the bias and scatter 
of picking errors by a series of GEV distributions whose parameters change slowly with decreasing SNR. In the next 
stage of this project, we will incorporate these non-Gaussian picking error models into our modified linearized 
location algorithm, hence utilizing a complete covariance matrix based on spatially correlated and SNR dependent 
errors. A hypothesis test will then be developed to verify when the linearized approach is valid. 
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